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I. SUPPLEMENTARY NOTES. CLONING RATE AND PROBABILITY 

DISTRIBUTION FUNCTION 

A. Evolution of walkers 

Hamilton's equations of motion with noise: q = p, p = — Vl^(q) + \/eri{t) can be written 
in a compact form: 

Xi = -Uij — + DiTJi (1) 

where x = (gi ■ ■ ■ qN,Pi ■ ■ -Pn), Di = (0 ■ • ■ 0, 1 ■ ■ ■ 1), cu = and repeted indices are 

summed over. The divergence of nearby trajectories is measured by the Lyapunov vector u 
which evolves as 

iti = -AijUj where Aij = ^ik (2) 
In terms of the unit vector v = j^, the norm of u is given by 

|u(t)| = \u{0)\e--^o^^^^m'^i (3) 

while V evolves as 

Vi = -AijVj + VkAkiViVi (4) 
The generating function of the largest Lyapunov exponent at time t consequently reads 

where the average is over trajectories given by ([T]) and (jlj), and over the noise realisations. 
Using standard methods, this average can be written 

(e"'*)traj. = y"dxidv, e-*(^--+"''^^-^'^)Po(xo,^o) (6) 

where Pq is the initial probability distribution and (xf, vj) is the end point of the trajectories. 
The operator Hpp reads 

e d"^ d dViq) d d 

2 dpi dqi dqi dpi dvi 

Equation ([7]) corresponds to an evolution: 

dP 

— = -{Hpp + av,AijVj)P (8) 
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It does not conserve probability: 

f dP f 

/ dxdv-^ = —a / dxdv ViAijVjP (9) 

This simply says that if one represents P{t) by a distribution of walkers, each one is cloned 
at every time step with a rate —aAijViVj, nothing but a times the stretching rate of |u(t)|. 

Clearly this derivation is not restricted to Hamiltonians of the form + V{q) and is 
valid also for the time dependent case. Furthermore, maps like the standard map can always 
be written as a continuous time dependent Hamiltonian evolution: 

kS 

q = p = — sin(27rg) for even time intervals 

(10) 

q = p p = for odd time intervals 

and the derivation above applies to them. 



B. Energy conserving noise 
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For a Hamiltonian of the form ^ + V{q), we implement an energy conserving stochastic 
dynamics by using a noise vector tangential to the energy surface: 

= E - ^) (11) 

where pj are independent white noises. When this equation is understood in the Stratonovich 
sense, it leads to microcanonical measure. 



C. Legendre transform 

Let us show at that time that by fixing a one can select - in the long time limit - 
trajectories with a desired Lyapunov exponent. By definition the generating function of the 
Lyapunov exponent is given by 

(e"'*)traj. = J dAp(A,t)e-"* (12) 

where p(A, t) is the probability density of A at time t. As is well-known, p(A) can be written 
for large times as a large deviation function: p(A, t) ~ e^^^^K By saddle point evaluation one 
then gets 

(e"^*)t.aj. ^ e*("^*+^(^*» (13) 
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where A* satisfies /'(A*) = —a. The dominating orbits indeed have A = A*, and (aA* + /(A*)) 
is the typical cloning rate. Thus, from A* and the cloning rate, which are given by the 
algorithm, one can reconstruct p(A). 
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' f-j '! The transition from order to chaos has been a major subject of research 

Q ■ since the work of Poincare, as it is relevant in areas ranging from the founda- 
tions of statistical physics to the stability of the solar system|l|, Along this 

2 ' transition, atypical structures like the first chaotic regions to appear, or the last 

^ regular islands to survive, play a crucial role in many physical situations. For in- 

^ I stance, resonances and separatrices determine the fate of planetary systems {^,2], 

' and localised objects like solitons and breathers provide mechanisms of energy 

^ , transport in nonlinear systems such as Bose-Einstein condensates and biological 



molecules js], Q]. Unfortunately, despite the fundamental progress made in the 

>• ■ last years, most of the numerical methods to locate these 'rare' trajectories are 
(N , 

! confined to low-dimensional or toy models, while the realms of statistical physics, 

^ ■ . . 

chemical reactions, or astronomy are still hard to reach. Here we implement an 

^ . eflflcient method that allows one to work in higher dimensions by selecting trajec- 
tories with unusual chaoticity. As an example, we study the Fermi-Pasta-Ulam 
nonlinear chain in equilibrium [61] and show that the algorithm rapidly singles out 
! the soliton solutions when searching for trajectories with low level of chaoticity, 
2 I and chaotic-breathers in the opposite situation. We expect the scheme to have 
natural applications in celestial mechanics and turbulence, where it can readily 
^ ' be combined with existing numerical methods. 

Complex dynamical systems are generically chaotic: two nearby trajectories diverge expo- 
nentially with time, at a rate given by the Lyapunov exponent: = |5x(0)| exp(tAorb.)- 
If chaoticity is sufficiently strong, almost every initial condition leads in the large-time limit 
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to the same exponent A = limAorb.(^)- Such an averaged description is not complete, as in 
many systems some trajectories will present atypically chaotic or regular behaviour extend- 
ing (in the future or in the past) for very long periods - even forever - depending on their 
initial conditions. 

A case in point is that of a planetary system : of the possible sets of masses and 

orbital elements compatible with observational error, only some correspond to a situation 
not involving a planet ejection in the recent past or near future. Another example, to 
which we shall return below, is that of nonlinear media, where spatially localised soliton and 
'breather' solutions (which have atypical Lyapunov exponents) are reached only from very 
specific initial conditions j^, 6|. 

In order to describe quantitatively the distribution of Lyapunov exponents of different 
trajectories, two closely related approaches, both inspired in thermodynamics, have been 
proposed in the past. In one approach, the sampling is over trajectories of fixed time and 
different initial conditions 8[, while in the other the dynamics is perturbed with a small 
random noise, and the sampling is made over different noise realisations 9|]. Lyapunov 
weighted dynamics (LWD) (a method proposed in the context of chemical reactions 10Ull|) 
can be modified to practically implement both these approaches - although here for brevity 
we discuss only the latter (see 'Method'). A population of walkers evolves in phase space with 
Hamiltonian dynamics, each one perturbed by an independent weak random force, which 
may be energy- conserving or not. The walkers carry with them a 'Lyapunov' vector that 
itself evolves as the separation between two nearby trajectories. After every time interval, 
walkers are cloned (or killed) with a rate proportional to {a times) the stretching rate of their 
Lyapunov vector. The end effect of the cloning is that orbits are now weighted according to 
their chaoticity. 

A positive (negative) value of a tends to favour orbits with large (small) Lyapunov expo- 
nent. In fact, this procedure counts the trajectories with a weight exp(aAorb.i) (See supple- 
mentary notes). Within the formalisms mentioned above, a and A are conjugate variables, 
just as inverse temperature and energy are in a thermodynamic problem: fixing a we obtain 
a typical value of A, together with its probability and the trajectories that contribute. 

In what follows we shall apply our method to several examples of increasing complexity. 
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Separatrices 

Separatrices are the phase-space frontiers between regions with different dynamical be- 
haviours. They play an important role because they are the cradle of chaos. In regular 
systems, walkers weighted with a = 1 can be shown to converge there, and to populate 
them uniformly. This is illustrated in figure 1 for the simple example of a double well with 
potential energy V{q) = + q^- 2g2. Particles can either librate within one of the wells, 
or oscillate between them. The separatrix is the curve emanating from the saddle point that 
limits these two regimes. Applying the dynamics with a = 1, the clones slowly diffuse in 
energy because of the noise until they converge to the separatrix, where they multiply more 
favourably. Clones that subsequently diffuse away die without leaving offspring. 




FIG. 1: Convergence towards the separatrix. 

The color code represents the value of the energy. Starting in one well (a, t ~ 3000), the walkers 
first diffuse in energy (b, t ~ 11450) until they reach the separatrix (c, t ~ 11725), where they 
settle (d, t > 12000). Two thousands walkers were run with e = 10^^ and a = 1. 

Let us now consider what happens as chaos sets in. We illustrate this with the classical 
example of the 'Standard Map': 

k5 

Pn+i =Pn- TT sin(27rg) = g„ + 5pn+i (1) 

/vr 

It represents the evolution of a free rotor regularly kicked with a constant force of fixed 
orientation j^. In the limit 5 — 0, it reduces to the evolution of a pendulum, and is thus 
integrable. For 5 > 0, it becomes chaotic, the more so the larger 6 and k. Applying the 
algorithm for a = 1 close to integrability (figure 2. a), the walkers concentrate on the unstable 



FIG. 2: The Standard Map 

The standard map trajectories are represented in orange. In blue, a thousand walkers at time 
T = 10,000 evolving with e = 10"^^. a) Quasi-integrable case: 5 = 0.41, k = 1. The map 
is very slightly chaotic, and the walkers converge for a = 1 to the homoclinic orbit emanating 
from the unstable fixed point. The inset is enlarged seventy-five times, b Secondary structures: 
(5 = 1, A; = 1. Several separatrices are revealed, starting from different regular islands with a = 1. 
c The last four regular islands in a strongly chaotic case [13] : = 1, k = 7.7a = —1. The 



insets are zoomed between 25 times (bottom right) and 150 times (top) and are centered around: 
(0.207; 0.09), (0.883; 0.09), (0.116; -0.09) and (0.8; -.09). 



manifold emanating from the fixed point, revealing the typical features of the homoclinic 
tangle. When chaos is increased — for larger 5 — many secondary structures become 
apparent (figure 2.b). Around the main resonant island, a stochastic layer is now the most 
chaotic structure of the phase space, and will consequently be the favoured target of the 
walkers. Nevertheless, starting from a given configuration, the walkers first converge to the 
closest stochastic layer, where they stay for some time. In other words, chaotic layers are 
the metastable states of the a = 1 dynamics. 

For a strongly chaotic system, we may wish to explore hidden regular structures 13] lost 
in the chaotic sea. Choosing a = — 1, we bias the measure in favour of regular orbits, thus 
revealing the — otherwise invisible — last resonant islands (see figure 2.c). 

Arnold Diffusion 

An integrable Hamiltonian system has a constant of motion for every degree of freedom. 
The phase space is filled with invariant tori on which the motion is quasi-periodic, with 
frequencies ui ■ ■ -un- Under a small perturbation, phase space is qualitatively unchanged, 
except for a chaotic region in the neighbourhood of 'resonances', defined as the points satis- 




FIG. 3: Arnold web and convergence. 

Arnold Web for the Hamiltonian (2) with fi = W~^, very close to integrability. Lighter regions 
have stronger chaoticity. In the inset: the evolution of a set of clones in three successive times, 
starting from the point indicated by the white cross. Yellow for < t < 160, 000, orange for 
160, 000 < t < 2, 000, 000 and red for 2, 100, 000 < t < 2, 200, 000. The large picture was composed 
by repeating the procedure shown in the inset for several initial conditions. Note that although 
the evolution of the clones resembles Arnold diffusion, it is in fact noise-driven, and hence orders 
of magnitude faster. 



fying the Diopha„t.„e Mation: E".-. ^ w.th „. integer In =ysteo,s with .note than two 

degrees of freedom, this 'chaotic web' leads to global diffusion in phase-space [IJ, 12|]. This 
phenomenon, predicted by Arnol'd in 1964 [l^, is called Arnol'd diffusion and the stochas- 
tic web is consequently called Arnol'd web. Localising such diffusive part of phase space 



is of great interest for many systems, 
physics or celestial mechanics. Laskar 



br instance in synchrotron experiments 



14| . plasma 



I5I developed a successful method, based on Fourier 
analysis, to map completely the regular part of phase-space. Here we are interested in going 
straight to the chaotic regions, which may be very thin and hard to find. Applying LWD 
with a > 0, the walkers first diffuse in phase space thanks to the stochastic noise, until they 
reach chaotic structures, where they multiply and settle. As a benchmark, we constructed 
the Arnol'd web for the Hamiltonian |l^: 
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FIG. 4: Finding resonances in higher dimensions. 

Main figure: a 6 degree of freedom system with Hamiltonian (2) and fi = 5.10"^, starting from a 
random initial condition. The two inflection points indicated by the arrows correspond to times at 
which the clones reached the resonances labeled (ni, 712, ns, n4, ns, ne, ^7). In the inset: the same 
plot for the dynamics in the inset of figure 3; yellow, orange and red correspond to the same three 
times shown there. 



N 2 

^ = E|+^-+i + ^^l^ / (2) 

~l ^ cos + cos t + A/ + 2 

Consider first the two-dimensional case N = 2. When /i = 0, this Hamiltonian is inte- 
grable, and the resonance lines are given by nipi + n2P2 = n^. For very small fi the chaotic 
set, shown in figure 3, is concentrated around the resonance lines. On the inset, one can see 
the path followed by the walkers that started in a random initial condition: they diffuse and 
settle in successive chaotic regions. The important point is that although the trajectory of 
the clones mimics Arnold diffusion, it happens in a timescale orders of magnitude shorter, 
as it is driven by noise. On the inset of figure 4 we show that each time the clones find a res- 
onance there is an inflection in the finite-time Lyapunov exponent. Figures like 3 have been 



obtained in 



16| by computing the finite-time Lyapunov exponents of trajectories starting in 



points on a grid in (^1,^2)- Here, instead, we run LWD as in the inset, starting from several 
initial conditions. 

In higher dimensions, the Arnold web is both difficult to represent and impossible to map 
exhaustively using a grid. There is no problem, however, in applying LWD as described 
above. To illustrate this, in figure 4 we show for the 6 degree of freedom version of the 
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FIG. 5: Finding solitons 

Energy-conserving (E = 1) dynamics of tiie Fermi-Pasta-Ulam chain with fixed ends, and strongly 
negative q, obtained starting from microcanonical equilibrium. The plot shows the position of the 
particles versus time. Several kinks bouncing from an end to the other of the chain are clearly 
visible. The Lyapunov exponent is half the typical value. (For clarity, the particles' average position 
has been arbitrarily displaced). 

running-time evaluation of the clones' Lyapunov exponent starting from a random config- 
uration. Using a separate programme, we have checked that the inflections indeed occured 
when the clones found the resonances indicated with arrows. Few-planet systems are clearly 
within reach, and LWD could be used to locate nearby chaotic trajectories in agreement 
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with experimental data, in the spirit of 
Breathers and Solitons 

As an example with many degrees of freedom, consider the Fermi-Pasta-Ulam chain, 
defined by the Hamiltonian: 

i ^ ' 

This chain is known to have localised soliton excitations, as well as the so-called 'chaotic 
breathers' 5|, la]. These excitations are unstable, but can be observed for some time starting 
from suitable initial conditions. Solitons configurations have atypically low chaoticity, while, 
perhaps more surprisingly, chaotic breathers are indeed very chaotic Q]. 

We have run an energy-conserving version of LWD on the chain starting from a micro- 
canonical equilibrium configuration. The results, for an energy density = 1, corresponding 
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Particle 



FIG. 6: Finding a chaotic breathers. 

Energy-conserving (energy density E = 1) dynamics of the Fermi-Pasta-Ulam chain with periodic 
boundary conditions and large, positive a, starting from microcanonical equihbrium. The configu- 
ration evolves towards a single chaotic breather, whose Lyapunov is three times the one of a typical 
equilibrium trajectory. The colour code corresponds to the energy of the site. 

to the 'transition region' j^, are shown in figures 5 and 6. When biasing to obtain a Lya- 
punov exponent three times larger than typical, we observe a chaotic-breather automatically 
emerging, thus showing that the large Lyapunov configurations are dominated by them. Con- 
versely, constraining the Lyapunov exponent to be a half of the typical value results in a 
few-kink configuration, as is clearly visible in figure 5. 
Conclusion 

Let us conclude by suggesting that the special trajectories described here play the same 
role as the reaction paths in physical chemistry, whose efficient detemination is currently 



a field of intense activity 



18| . In all cases one has to deal with events that are important 



but happen rarely, so they are best studied with a dynamics that is biased to enhance 
the probability of their occurrence - although of course one is ultimately interested in the 
properties of the true, physical dynamics. 

Method A population of walkers with phase-space positions = (q*^, p*^) is evolved 
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with the standard Hamiltonian dynamics q** = p*^, p** = — Vl^(q'*) + ^/er]{t), where is a 
white noise of unit variance. Every clone has a 'companion' which starts a small distance 
5x^(0) away, and evolves with the same noise. After a time interval 6t, when the positions of 
clone and companion are [{x."-{6t),x.^{6t) + 6x."-{6t)], and the separation ratio Pa = [ \\s^a(o)\\ ) 



i) their mutual separation is normalised — > 6x."-/pa, ii) if (Pa)° > 1 the clone a is copied 
(i.e. replaced by two clones) with probability (pa)" — 1, while if (pa)" < 1 it is killed with 
probability 1 — (pa)"- Each clone-companion pair subsquently evolves with a different noise 
realisation. An overall cloning rate is then applied to keep the population constant. In the 
energy-conserving variant, which we applied to the particle chain, the only difference is that 
the noise vector is tangential to the energy surface. 
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